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In this paper we investigate how the halo mass function evolves with redshift, based on 
a suite of very large (with Np = 3072'^ — 6000'^ particles) cosmological N-body simula- 
tions. Our halo catalogue data spans a redshift range of z = — 30, allowing us to probe 
the mass function from the dark ages to the present. We utilise both the Friends-of-Friends 
(FOF) and Spherical Overdensity (SO) halofinding methods to directly compare the mass 
function derived using these commonly used halo definitions. The mass function from SO 
haloes exhibits a clear evolution with redshift, especially during the recent era of dark energy 
dominance (z < 1). We provide a redshift-parameterised fit for the SO mass function valid for 
the entire redshift range to within ~ 20% as well as a scheme to calculate the mass function 
for haloes with arbitrary overdensities. The FOF mass function displays a weaker evolution 
with redshift. We provide a 'universal' fit for the FOF mass function, fitted to data across the 
entire redshift range simultaneously, and observe redshift evolution in our data versus this fit. 
The relative evolution of the mass functions derived via the two methods is compared and we 
find that the mass functions most closely match at z = 0. The disparity at z = between 
the FOF and SO mass functions resides in their high mass tails where the collapsed fraction 
of mass in SO haloes is ^ 80% of that in FOF haloes. This difference grows with redshift so 
that, by z > 20, the SO algorithm finds a ^ 50 — 80% lower collapsed fraction in high mass 
haloes than does the FOF algorithm, due in part to the significant over-linking effects known 
to affect the FOF method. 

Key words: Large-scale structure of the universe: cosmology:theory — dark matter — 
galaxies:haloes — galaxies:high-redshift — methods: numerical 
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1 INTRODUCTION 

The abundance of dark matter haloes - the halo mass function 
- plays an important role in cosmology due to its sensitivity 
to a number of important parameters including the matter den- 
sity of the universe, Q.m, the Hubble parameter, h, the spectral 
index of the primordial power spectrum, Us, and the dark en- 
ergy equation of st ate (e.g. lHolder et alJboOltlHaiman et al I2OOII : 
IWelIeretalJl2002l) . Whilst the mass function is curre ntly diffi- 
cult t o pin down with much precision observationally ( Eke et all 
l2006l) . it is relatively easy to pr obe via N-body cosmological simu- 
lations (see lLuki"c eTaD 12007!) for a review of older work, more 
recent studie s includ e Tinker etal. (20oi); Crocce et alj 1 2O10h: 
ICourtin et all ( l201lh : jBhattacharva et alj i201lh ; lAnguIo et al.l 
( l2012h ) and can, to a n extent, be understood an d modelled through 
analytic arguments jPress & Schechterl 1 1974 ISond et al.1 Il99ll : 
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Lee & Shandarinlll998l:ISheth & Tormeij|2002l; 


Maeaiore & Riottd 


2010l;ICorasaniti & Achitouvll201 ll; Lim & Lee 


2012h. 



The halo mass function is employed widely in astronomy: at 
low redshifts it is used in statistica l analyses of cluster su rveys to 
constrain cosmological parameters ( ISmith&Mariarj201 l|) as these 
large objects probe the high mass end of the function (for a re- 
cent review paper on the formati on of galaxy clusters , inclu ding 
the role of the mass function, see iKravtsov & Borganil j2012l) ). It 
is a key component in studies employing correlation functions be- 
tween galaxies, as the halo-halo term directly depends on the mass 
function (see for example Coorav & Sheth (2002) and references 
therein). At higher redshifts it is applied in modelling the process 
of reionisation, which proceeded between z ~ 6 — 20 from sources 
residing in dark matter haloes including population-Ill stars, early 
galaxies and accreting black holes. Any significant observed devi- 
ation from the mass function predicted by the ACDM cosmolog- 
ical model would create tension in our current understanding of 
structure formation in the universe. For example we are able to put 
bounds on how massive the largest observable clusters in the visible 
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unive rse should be from the mass function (see [Harrison & ColesI 
and references therein), the discovery of larger clusters would 
require explanation. 

Much has been been written on the topic of the mass func- 
tion by groups working with N-body simulations. To date this has 
largely focussed on low redshift {z < 2) studies - due to the dif- 
ficulty in simulating haloes at high redshifts with any great accu- 
racy - a nd haloes derived using the Friends-of-Friends (FOF) al- 
gorith m (P avis et aL 1985). There are notable exceptions, for ex- 
ample |R^d_etaL. ( 2007.) focuses on the FO F mass function at high 
redshifts (z > 10) and [Tinker et al.l ( 1200 Sh give a precision mass 
function f or haloes derived using the Spherical Overdensity (SO) 
algorithm dLacev & Cole 1994). There is currently no robust mass 
function available for haloes derived using the SO algorithm at red- 
shifts higher than z ~ 2 which represents a significant gap in the 
literature that we address in this work. 

The majority of mass function studies performed to date 
in simulations have employed the FOF algorithm, not least be- 
cause it is easier to implement and often computationally ch eaper 
than the SO algori t hm. D espite this, [Tinker et all ( |2008|) and 
iKravtsov & Borganil ( 1201 2h argue that an SO halo definition is 
much more appropriate for correlating halo masses to X-ray, 
Sunyaev-Zeldovich (SZ) effect and optical observables. One of the 
purposes of this paper is to provide a useful SO mass function, cali- 
brated to a ACDM cosmology, that covers as broad a redshift range 
as possible. We also calculate a universal FOF mass function and 
investigate how our data evolves with redshift against it. 

The question of whether the mass function is 'universal' - in- 
dependent of redshift or cosmology, when expressed in suitable 
units - has been answered recently in the negati ve by various 
authors ([ke ed et al. 2007; Tinker et al. 2008; Croc ce et alj|201(]| ; 
ICourtin et al.ll201 ll ; iBhattacharva et alj|201 ih despite there being 
studies prior to these that found little evidence for d epartures from 
universality djenkins et alj|2()0ll; iLukic et alj|2007l) . It was origi- 
nally noted bv I Jenkins et al. I ( I2OOII)" that in order to produce a uni- 
versal mass function from FOF haloes one needs to take a constant 
linking length - the key parameter that the FOF algorithm uses, see 
§ [3] - across all redshifts or cosmologies being considered. How 
appropriate this is is questionable as it has been observed that an 
important factor in creating FOF-based halo catalogues is the sup- 
pression of a systematic effect inherent to the algorithm known as 
'overlinking' which occurs when two haloes are linked together 
by a bridge of particles. When this occurs the resultant objects are 
in no way physically interpretable as virialised collapsed density 
peaks and yet are almost always unavoidably includ ed in FOF mass 
funct i ons. Previous author s tPavis et al.it 1985. ; .Ma & BertschingeJ 
ll995l;ICole & Lacevlll99^ ; lLukic et al.ll2009h have calculated that 
15-20% of all FOF haloes calculated with the standard linking 
length of 0.2 are objects that have been bridged together in this 
manner and that larger linking lengths lead to more overlinking. It 
is not clear to what extent overlinking impacts on the universality of 
the FOF mass function. Is the function closer to universal because 
the FOF haloes are overlinked in ~ 20% of cases, or does the over- 
linking have a nega tive i mpact on the univ ersal result? Studies by 
ICourtin et al] ( l20Ilh and lMore etal] ( I2OIII) suggest that a suitable 
choice for the linking length parameter leads to a universal mass 
function and that the effect of overlinking is to distort this result. 

[Tinker et alj ( l2008h observed a much clearer departure from 
universality for S O haloes than for FOF haloes, a fact that 
iMore et al] ( l20Ilh propose is due to the effect of taking a fixed 
overdensity criterion - the key parameter of the SO algorithm, see 
§ [3]- rather than specific overdensity criteria bespoke to a given 



redshift or cosmology. One is therefore left with a choice: either 
adjust one's halo definitions to attempt to produce a universal mass 
function or keep one's halo definitions fixed and expect the func- 
tion to deviate from a universal form. In this work we choose the 
latter. Deciding upon a suitable 'universal' redshift-dependent halo 
overde nsity or linking length para meter is stil l an ope n question (al- 
though ICourtin et al] ( I2OIII) and iMore et alj ( I2OI II) go some way 
towards providing a solution), due to complicating factors such as 
the triaxiality of haloes in the SO case - a purely spherical halo is 
an ideal that is never realised - and overlinking in the FOF case. 
We shy away from demanding a universal mass function and pin 
our halo definitions down firmly in § [3] 

This paper is laid out as follows. In §[2]we describe the simula- 
tions we have used to construct the mass function. In §|3]we outline 
the halo-finding methods we employed and the definitions of what 
we are referring to by the term 'halo'. In §|4]we present our mass 
function results and fits. In §[5] we present a summary of the work 
and a discussion of various points that arise. 



2 SIMULATIONS 

The simulations undertaken in this work are summarised in Ta- 
ble □ They follow between 3072^ (29 billion) to 6000^ (216 bil- 
lion) particles in a wide range of box sizes from 11.4 /i^^Mpc 
up to 6 ft^^Gpc. Spatial resolutions range from 0.18 /i^^kpc to 
50 /i~^kpc while particle masses range from 3.6 x lO^/i"^M0 to 
7.5 X lO^''/i~^M0. This allows dark matter haloes to be resolved 
in a very large mass range, with a low end of 7.3 x lO*/i^^M0 
(for a 20 particle halo) in our smallest volumes, and no upper 
limit in the halo mass, since our largest volumes approach the 
size of the observable universe. All simulations were performe d 
using the CubeP^M N-body code dHamois-Deraps et al] l2012h . 
The CubeP^M code calculates short-range direct particle-particle 
forces combined with a long-range PM force calculation making 
it a P'^M (particle-particle-particle-mesh) code. It is massively par- 
allel and runs efficiently on either distributed- or shared-memory 
machines. This is achieved via a cubical, equal-volume domain de- 
composition combined with a hybrid OpenMP and MPI approach. 
CubeP'^M scales well up to thousands of pr ocessors and to date has 
been run on up to 21,976 co mputing cores dlliev et al.ll2008ll20I(]|; 
iHarnois-Deraps et alj|201 j) . We note that such large simulations 
yield very large amounts of data, with just the particle data amount- 
ing to between 700 GB (3072^) and 4.7 TB (6000^) per time-slice, 
providing significant challenges in the data handling and analysis, 
as we discuss below. 



2.1 Cosmology 

We base our simulations on the 5 -year WMAP results 
dDunklev et al]|2009l ; IXomatsu et alj|2009h . With the exception of 
two of our runs (the 1 /i~^Gpc and 3.2 /i~^Gpc boxes) the cosmol- 
ogy used for t h e sim ulations was the 'Union' combination from 
iKomatsu et al] d2009l) . based on results from WMAP, baryonic 
acoustic oscillations and high-redshift supemovae; i.e. fim = 0.27, 
Q.A = 0.73, h = 0.7, = 0.044, gs = 0.8, Us = 0.96. The 1 
/i^^Gpc and 3.2 ft^^Gpc boxes were based on the slightly different 
'Alternative' combination from the same paper; i.e. Q,m = 0.279; 
f^A = 0.721, h = 0.701, ^6 = 0.046, erg = 0.817, Us = 0.96. 
The power spectrum and transfer funct ion used for settin g initial 
conditions was generated using CAMB dLewis et al.ll2000l) . 
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Table 1. N-body simulation parameters. Background cosmology is based on the WMAP 5-year results. Minimum halo mass and the redshift of formation of 
the first halo are base d on a 20 particle minimum . Resolution here refers specifically to the softening length, above which the interparticle force is the exact 
Newtonian force fsee lHamois- Derapsetal.l)2012l) for details). See text for details regarding the quantities 'P{k)max and 5™'' /Ap. 
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2.2 Initial Redshift and Initial Conditions 



3 HALO FINDING 



The CubeP M code uses first-order Lagrangia n perturbation the - 
ory (ILPT), i.e. the Zel'dovich approximation i lZerdovicd[l970h . 
to place particles in their initial positions. The initial redshift when 
this step takes place is given for each simulation in table 1. 

There is debate regarding the suitab ility of the Zel'dovich ap- 
proximation in setting initial conditions I Crocce et al. 2006, 2010; 
Lukic et al.l2007l : lTinker et alj2008l:lKnebe et al.ll2009l : IReed et al.l 



2012h . Care must be taken when setting a choice for the ini- 
tial redshift so as to ensure that artefac t s (tran sients) are min- 
imised. A detailed study by iLukic et ^ ( l2007h emphasises the 
need for a suitable choice for the initial redshift as applied to 
a particle- mesh (PM) code (MC - see iHeitmann et al. I l2005l : 
iHeitmaniT et al. 2008) using ILPT. iReed et al.l ( 1201 2h recommend 
the use of second-order Lagrangian perturbation theory (2LPT) in 
setting initial con ditions b ased on their study using two tree-based 
codes: Gadget-2 jSpringel 2005) and PKDGRAV ( Stadel 2001). 



iLukic et^ ( l2007h proposed two criteria to guide the choice 
of initial redshift for simulations employing ILPT: (1) whether 
the amplitude of the initial power spectrum modes in the box are 
in the linear regime, and (2) whether the initial particle displace- 
ment is small enough so that the particle grid distortion is relatively 
small. The first criterion sets a minimum for the initial redshift 
based on the requirement that the dimensionless power spectrum, 
Vik) = it3p(fc)/27rVbeJess than some arbitrary, small value at 
the initial redshift ( iLukic et alj ( 12007 ) used 'P{k)max < 0.01). As 
smaller boxes for a given particle number probe higher values of 
fc, earlier initial redshifts are required for smaller boxes. Similarly, 
for a fixed box size, simulations with greater numbers of particles 
also require earlier initial redshifts. The second criterion places a 
m ore stringent bound on the initial redshift. The example discussed 
m iLukic et alj ( |2007|) was that the particles were displaced, on av- 
erage, by no more than S^^" — 0.3Ap, where Ap = Lbox/np 
is the interparticle spacing and 5™" is the root-mean-squared dis- 
placement of the particles. In Table 1 the values of 'P{k)max and 
S™" / Ap for our simulations are given. We see that for all of our 
simulations both of these criteria are satisfied. 



3.1 Motivation 

In the hierarchical picture of structure formation peaks in the linear 
density field of the early universe grow into highly non-linear struc- 
tures via gravitational attraction. A popular formalism for describ- 
ing this process is the identification of the non-linear, high den- 
sity peaks as dark matter haloes. These haloes, appropriately de- 
fined, are concomitant with virialised regions in the universe where 
galaxies reside. The hierarchical growth of structure in this picture 
can then be viewed as a series of halo mergers, commencing with 
very small haloes at early times, and continuing all the way to the 
present, with haloes growing larger via both mergers and the ac- 
cretion of smaller haloes. In cosmological simulations it is usual to 
attempt to track the progress of the growth of structure by identify- 
ing haloes within the simulation volume across a range of redshifts. 



3.2 Halo Finding Approaches 

There are two main approaches to finding haloes in cosmological 
simulations. The first is the spherical overdensity (SO) method. In 
this method haloes are assumed to be spherical. The extent of a 
halo is governed by a free parameter, the overdensity criterion. A, 
which is a cut-off in density with respect to some background den- 
sity (typically either the matter background density or the critical 
density of the universe). In the SO algorithm spheres are grown 
from a central location until the enclosed overdensity of the sphere 
is equal to the overdensity criterion. The main differences between 
halofinding codes that use this approach lie in how the centres of 
the candidate haloes are identified. 

The second approach is the Friends-of-Friends (FOF) algo- 
rithm. This algorithm is also based on one parameter, b, the linking 
length parameter. The algorithm finds haloes that contain particles 
that are within bAp of at least one other particle in the halo. While 
the SO algorithm produces a halo that is by definition spherical, 
the FOF algorithm creates haloes that are arbitrarily shaped. For 
a recent comparison project on halofinding codes see lKnebe et al.l 
(2011 ). 

We employ three halo-finding codes in our analysis: 
CubeP^M's own on-the-fly SO halofinder (hereafter 'CPMSO') 
dHamois-Deraps et alj I2OI2I) , the Amiga Halo Finder (hereafter 
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'AHF') join et alJl2004lKnollmann & Knebdlioosl) . and the FOF 
halofinder from the Gadget-3 N-Body cosm ological code ( an up- 
date to the pubhcly available Gadget-2 code l lSpringei2005h ). 

The CPMSO halofinder utilises a fine mesh from the 
CubeP^M code (a mesh with spacing of Ap/2) to identify local 
peaks in the density field. The code first builds the fine-mesh den- 
sity using either Cloud-In-Cell (CIC) or Nearest-Grid-Point (NGP) 
interpolation. It then proceeds to search for and record all local den- 
sity maxima above a certain threshold (typically set to 100 above 
the mean density) within the physical volume. It then uses quadratic 
interpolation on the density field to determine more precisely the 
location of the maximum within the densest cell. The halo cen- 
tre determined this way agrees closely with the centre-of-mass of 
the halo particles. Each of the halo candidates is inspected inde- 
pendently, starting with the highest peak. The grid mass is accu- 
mulated in spherical shells of fine grid cells surrounding the max- 
imum until the mean overdensity within the halo drops below A. 
While the mass is accumulated it is removed from the mesh, so 
that no mass element is double-counted. This method is thus inap- 
propriate for finding sub-haloes as within this framework they are 
naturally incorporated in their host haloes. Because the haloes are 
found on a grid of finite-sized cells and spherical shells constructed 
from them, it is possible, especially for the low-mass haloes, to 
overshoot the target overdensity. When this occurs we use an an- 
alytical halo density profile to correct the halo mass and radius to 
the values corresponding to the target overdensity. This analytical 
density profile is given by the Trunc ated Isothermal Sphere (TIS) 
profile ( Sh apiro et al.ll 199^ ; llliev & Sh apiro 2001) for overdensi- 
ties below ~ 130, and 1/r^ for lower overdensities. The TIS den- 
sity profile has a similar outer slope (the rele vant one here) t o the 
Navarro, Frenk and White (NFW) profile ( Nav arro et al .Il997h .but 
extends to lower overdensities and matches well the virialization 
shoc k position given by the Bertschinger self-similar collapse solu- 
tio n teertschingeil 198a). For fur ther details on the CPMSO method 
see lHamois-Deraps et al.l ( [20121) . 

Details for AHF can be found in iKnollmann & Knebd |2009|) . 
The algorithm identifies density peaks using a nested set of grids 
that are finer-grained in regions of higher density. Haloes are then 
identified by collecting particles together that are contained within 
isodensity contours on each level of the grid hierarchy. This allows 
AHF to identify sub-haloes within host haloes as the algorithm pro- 
gresses from high resolution grids to lower resolution ones. AHF 
then removes particles that are unbound and recalculates halo prop- 
erties based on the remaining bound particles. 

The specifics of the FOF halofinder packaged in with the 
Gadget-3 code currently have not been detailed in any pu blica- 
tion but the algorithm itself is outlined in iDavis et all ( 1 19851) . The 
main difference in the algorithm that exists in the Gadget-3 version 
is that the code is parallelised for distributed-memory machines. 
Specifically, haloes are found in local subvolumes of the simula- 
tion assigned to individual MPI tasks (created using the Gadget-3 
domain decomposition which utilises a spac e-filling P eano-Hilbert 
curve - for details see the Gadget-2 paper 'Springel' and 
then haloes that extend spatially beyond the edges of the subvol- 
umes are linked together in a final MPI communication step. We 
have altered the Gadget-3 code to read CubeP'^M's particle output 
format and significantly reduced its memory footprint by stripping 
away extraneous data structures. 

Finally, we note that due to limitations in the scaling of the 
codes with processor numbers, large memory footprint, and incom- 
patible data structures, in order to apply the AHF and FOF algo- 
rithms to our data it was necessary to split the simulation time- 



slices into a number of subvolumes and run the halofinding al- 
gorithms on each subvolume independently. Each subvolume in- 
cluded a buffer zone which overlapped with the neighbouring ones, 
for correct handling of haloes straddling two or more sub-regions. 
We then stitched the subvolumes back together to create the final 
AHF and FOF halo catalogues, removing any duplicated structures 
in the overlapping buffers. Although somewhat more expensive 
than applying each halo finder directly on the full data, this ap- 
proach allows the handling of much larger amounts of data than 
otherwise possible and provides additional flexibility in terms of 
computational resources needed for post-processing. 



3.3 Halo Definitions 

In this work we choose the overdensity criteria for identification of 
spherical overdensity haloes to be Ai78, i.e. an overdensity of 178 
times the background matter density. This is a common choice mo- 
tivated from the top-hat model of non-linear c ollapse in an Einstein 
de-Sitter (EdS) universe ( lGunn&Gottll972l) . The overdensity cri- 
terion is usually taken to be an overdensity with respect to either 
the background matter density or the background critical density. 
For an EdS universe, i.e. a universe containing only coUisionless 
matter, the two are the same. For a ACDM universe, at late times, 
the two deviate from each other due to the increasing dark energy 
component that contributes to the critical density of the universe but 
not to the matter density. Care must be taken to compare like-for- 
like overdensities, especially during the A-dominated epoch (note 
that all references in this work to overdensities refer to those with 
respect to the background matter density). Any given non-linear 
overdensity criterion can be mapped onto a threshold value in the 
linear regime. For Airg this corresponds to a linear overdensity of 
5c — 1.686 in an EdS universe. In a ACDM universe the value 
of 5c evolves to 1.674 al z = due to the influence of dark en- 
ergy. In the case of the FOF halofinder we follow various previous 
authors Jjenkins et al.|[200ll: iReed et alj|2003l |2007| ; ICrocce et all 
l20ld : ICourtin et al.l 



2011 



Angulo et al 



20121) and use a linking 



length of 0.2 for our analysis. 



3.4 Halo Mass Redsliift Dependence 

There are a number of factors that influence the growth of halo 
mass. They can be divided into 'physical' mass growth, i.e. the 
growth of the haloe s via mergers and accretion, and, following 
the nomenclature of Diemer et al] j20I2l) . 'pseudo' mass growth. 
For a recent study on the relative effects of physical and pseudo 
mass growth on the SO mass function between 2 = — 1 see 
iDiemer et alj ( 1201 2h and references therein. An important point 
from this study is that a significant fraction of SO halo mass growth 
between z = — 1 is due to pseudo-evolution of halo mass. This is 
growth of halo mass, defined in the SO sense, due to the evolution 
of the background density of the universe. As the background mat- 
ter density of the universe decreases with the universe's expansion, 
the radius of a given SO halo will grow as it requires a lesser en- 
closed physical density to meet the overdensity criterion of A times 
the background density. In the FOF case there also exists a pseudo 
mass growth: as the linking length is defined as a constant length in 
co-moving coordinates its physical length increases as z decreases. 
The implication of this is that for a halo wit h a static profile the ex- 
tent of an FOF halo will increase with time. lFakhouri&Ma|j20I0l) 
split FOF halo growth into an accretion component and a diffuse 
component, the latter containing an element of pseudo growth. Cur- 
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rently there is no equivalent study to the lOiemer et alj ( 1201 2|) work 
on pseudo mass growth in SO haloes for FOF haloes. 

In addition to the pseudo evolution of halo mass there is the 
evolution of 5c (and A) due to dark energy dominance at low red- 
shifts. This has well-studied implications for the SO halo mass defi- 
nition jLahav et alJl99ll : lLacev & ColJl993l ; lEke et aljl996h . Be- 
cause of this 5c(z) is often us ed as a parameter in mass functions 
1974; Shet h & Tormerjl 20021 : iReed et al.ll2003l 
uOUi) . In this work we provide an SO mass 



I Press & Schechtei 



l2007l : ICourtin et al 



function fit for haloes in an ACDM universe that is based on Ai7g 
and do not include 5c.{z) as a variable in the fit. This is motivated 
by the fact that we parameterise redshift evolution occurring due 
to dark energy, pseudo mass growth, and systematic effects such 
as halo mergers (which occur at different rates in different epochs) 
and non- sphericity of haloes. These latter effects are not captured 
by the evolution of 5c.{z). We also note that a baryonic c omponent 
will h ave an effect on the mass function (estimated bv ICui et al.l 
( l2012h to increase the mass function by < 3% for a A200 halo 
definition) not captured in this work. 

An issue that pseudo mass evolution raises is that any redshift 
evolution of the mass function (whether expressed in terms of mass 
or Incr^^, see §0 is due to an interplay between halo definition 
and the physical increase of halo mass with time. In order to link 
halo definitions to observations and compare our results to previ- 
ous studies on the halo mass function we adopt the common halo 
definitions above and attempt to capture redshift evolution based 
on them. 



4 THE HALO MASS FUNCTION 

The mass function can be expressed in a number of ways. Here we 
use the multiplicity function, /(cr, z), which represents the frac- 
tion of mass that has collapsed to form haloes per unit interval in 
lno"~^. This definition uses cr"^ as a proxy for mass, where is 
the variance of the linear density field, given by: 



cr^(M,z) = 



27r2 



k^P{k)W^{k-M)Ak 



(1) 



where P{k) is the power spectrum of the linear density field, 
W{k\ M) is the Fourier-space representation of a real-space top- 
hat filter containing mass M - assuming that the top-hat sphere en- 
closes a region that contains a mean density identical to that of the 
univers e - and D(z) i s the growth factor, normalised to the unity at 
z = Q jPeebleslll993l) . The radius in real-space of the filter can be 
set using an overdensity criterion. A, as per the SO algorithm. The 
relationship between a and halo mass for our two WMAP-5-based 

cosmologies is shown in Figure[T] 

W e define the halo multiplicity function as jjenkins et al.l 
l200lh : 



M dN{M,z) 



(2) 



where N{M, z) is the mass function proper, that is to say the num- 
ber count of haloes with a mass less than M per unit volume, and 
Pm is the mean matter density of the universe. 

4.1 Mass Binning 

The halofinders calculate values for the masses of all resolved 
haloes. It is then necessary, in order to construct /(a, z), to count 
the number of haloes in mass bins before converting the masses 
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Figure 1. The relationship between the vaiiance of the Hnear density field, 
cr, and halo mass fo r two ACDM c osmologies. 'Cosmology A' refers to 
the 'Union' model of lKomatsu et al.l 12009.) with flm = 0.27, Ua = 0.73, 
h = 0.7, = 0.044, as = 0.8, ris = 0.96. 'Cosmology B' refers 
to the 'Alternative' model of lKomatsu et alj j2009l) with Q^n = 0.279; 
Ua = 0.721, h = 0.701, = 0.046, ag = 0.817, = 0.96. 



into values of Incr 
catalogues we have: 



Rather than equation 2, for simulated halo 



pm{z) AM diner- 1 



(3) 



where care must be taken in choosing the width of the bins in mass 
as this can potentially result in a so urce of systematic error. For a 
detailed evaluation of this error see iLukic et al. I j2007l) , who note 
that the error is not significant so long as the bin widths, AlogAf, 
do not exceed 0.5. We use bin widths for our analysis that remain 
constant in logAf with AlogA/=0.I6. and assign the mass of a 
given bin to the average of all the haloes in the bin rather than the 
bin centre. 



4.2 Error Treatment 

Throughout this paper we use l-a Poisson error bars as defined in 
lHeinrichl ( [2003l) : 



a± = \ N haloes + 



4 2 



(4) 



The main advantages of using error bars in this form are: (I) for 
small numbers of haloes the asymmetry of the bars reflects the 
asymmetry of the Poisson distribution; (2) for a halo count of 1 
the lower edge of the error bar does not reach zero. Note that for 
large numbers of haloes the errors tend to the expected \/N form. 
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4.3 Data Treatment 

In order to combine data across our simulations certain system- 
atic effects need to be accounted for. These include: accounting for 
the finite volumes of our simulations; adjusting the FOF haloes for 
a systematic overestimation in mass for haloes sampled with low 
particle numbers; the question of whether to remove or include sub- 
haloes from our AHF data; and the choice of a lower limit to the 
number of particles a halo contains. The mass functions presented 
here are (1) corrected for finite volumes for simulations with a box 
size of less than 425 /i~^Mpc, (2) adjusted, in the case of our FOF 
haloes, to account for mass ov erestimation using the empirical cor- 
rection of IWarren et alj ( l2006h . (3) constructed based only on host 
haloes (in the case of data from the AHF halofinder the host haloes 
have masses that include any sub-haloes they contain) and, (4) have 
a minimum of 1000 particles in each halo. We now discuss each of 
these points in detail. 



4.3.1 Finite Volume Correction 

The volume of space modelled by a cosmological simulation is al- 
ways finite. However when we speak of the mass function we typ- 
ically refer to a 'global' mass function, i.e. one that would corre- 
spond to an infinitely large volume. There is therefore a disconnect 
that needs to b e bridged between simulations and an ideal, global, 
mass function ( ISirkoll2005l ; [Bagla & Prasa3l2006l : |Power & Knebd 
l2006h . There exist a number of ways to address th is issue. In 
this w ork w e adopt an approach recently employed bv lLukic et al.l 
( l2007h and iBhattacharva et all j201 lb . We proceed by assuming 
that mass function universality holds strictly in the sense that the 
functional form is the same for both global and local volumes. This 
is similar to assuming that the mass function is universal across 
cosmologies that differ in their mass-a relations, as the effect of a 
finite box size is to set to zero the amplitude of any density fluctu- 
ations on a scale greater than the box size. There is evidence that 
the mass function is not universal (see §[5]l; that its functional form 
has a weak dependence on both redshift and cosmology. Therefore 
any finite volume correction we make based on the assumption of 
universality is an approximation. We adopt the approach here de- 
spite this, as it is relatively straight-forward to apply and it brings 
our data into better agreement across difference box sizes. We also 
briefly discuss below other possible methods that could be applied. 

The approach is as follows. Assuming equation |2]refers to the 
mass function in an infinite simulation volume we can re-write it 



10 



dM 

and now for a finite simulation volume we can write 

dN' " 
dM' 



Pm ,, > diner' 



(5) 



(6) 



dlogM' 

where now a'{M') is determined by the discrete power spectrum of 
the simulation in question. Also M' is a function of AI defined such 
that a{M) = a' {M' (M)). The assumption that the mass function 
is universal allows us to say that the multiplicity functions, /, in 
equations |5] and |6] are identical. This then leads to the relation: 

,dM'(M) 



dN = dN 



dM 



(7) 



We now require some method of connecting a'(M') and a( M). 
The extended Press-Schechter formalism aBond et al.l Il99ll) ap- 
proximately connects a'{M') and a{M) via: 
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Figure 2. The finite volume correction for small boxes. The global M ~ a 
relation is adjusted to give a relation based on the volume of the simula- 
tion box in question via equation [8] The (approximate) assumption of mass 
function universality is then employed to adjust the masses and numbers of 
dark matter haloes via equations|2]and[8] 



(8) 



where cr 



R{box) 



is the variance of fluctuations in spheres that contain 
the simulation volume. To align this to the simulations we make an- 
other approximation and equate the spherical volume to the cubical 
simulation volume. The resulting mass-a relations at 2 = for the 
1 14 /i^^Mpc and 20 /i~^Mpc boxes are shown in Figure[2] 

The steps in the volume correction are then as follows: (1) 
Calculate cr^^j^^) for the simulation box in question. (2) For each 
bin in a calculate an adjusted bin a value, (Jcor, using a^^^ = — 
'^R(box)- (3) Adjust each halo count using equation [7] where M is 
the mass that corresponds to a and M' the mass that corresponds 

to acor- 

For comparative purposes we now discuss other so lutions that 
have b een employed to solve the finite volume issue. iReed et al] 
1 I2OO7!) performed a number of different N-body simulations and 
calculated the range of a for each one using the input power spectra 
of the simulations and the relationship between a and mass for a 
finite box in the discrete case: 



D\z)J2\S^\'w\k;M) 



(9) 



where |5k| is the linear amplitude of the Fourier modes in the simu- 
lation at 2 = 0. This approach has the advantage of both correcting 
for finite volumes and also compensating for cosmic variance. One 
drawback to it though is that each realisation has a different mass-a 
relationship. We did not adopt this approach as it is a method more 
suited to multiple simulation runs, where the mass-a relations can 
be averaged over to produce a mass function fit. We did not have 
the luxu ry of repeating our runs due to the large sizes of the sim- 
ulations. iBagla et all ( I2OO9I) choose to make no specific correction 
(although they remove any data points that are affected by by more 
than a threshold level of error - 10% in the number counts - from 
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their analysis in a similar manner to [Tinker et alj ^20081) '). This de- 
cision was motivated by the observation that the mass function is an 
unknown function that is deduced from simulation data so to place 
any a priori constraints upon it - such as a universal functional 
form for global and local mass-cr relations - is undes irable. This 
approa ch was applicable in the study undertaken by iBagla et al] 
l l2009l) as the box sizes used for their simulations were all relatively 
large 256 /i^^Mpc), for our smaller boxes not correcting for fi- 
nite volumes would lea d to undesirab l e syste matic disco ntinuities 
across our simulations. lYoshida et alj ( l2003h and iBagla & Prasad 
( l2006h replaced equation 1 with: 



27r2 



k'^P{k)W^{k;M)dk 



(10) 



/27r/L 

which take s a cu t-off in low k modes at the size of the box. 
iLukic et al.l bOOTh note that a correction of this form has a depen- 
dence on the accuracy of the mass functions at redshifts greater 
than z = 5. They also note that in comparison to their correction 
method it exhi bits offsets and shape c hanges across different box 
sizes. Finallv. iBarkana & Loebl j2004h utilised the extended Press- 
Schechter formalism to create a volume adjustment by taking the 
Sheth-Tormen mass function for a global volume and applying a 
correction to it to create a mass function suitable for a smaller vol- 
ume. As we are not adopting a mass function in the Sheth-Tormen 
form this approach was inappropriate here. 

4.3.2 Warren Correction to FOF Haloes 



IWarren et all ( |2006|) proposed a simple correction to the masses of 
FOF haloes that alleviates a systematic error in halo masses cal- 
culated via the FOF algorithm at low particle counts. This correc- 
tion was devised bas ed on analy s is of F OF haloes at z = 0, but 
has been checked bv lLukic et al] ( l2007t) for FOF haloes at higher 
redshifts. The FOF algorithm overestimates the masses of haloes 
when there are low numbers of particles sampling the haloes. 
IWarren et al.l (j2006) proposed the following correction to particle 
counts in haloes: 



N,orrected = N{l-N~°''') 



(11) 



We have adopted this correction in this work and all FOF data 
shown includes it. As we have also adopted a cut-off in particle 
number of 1000 the maximum effect this correction has on the 
masses of our FOF haloes is ~ 2%. 

There is deba t e rega rding the appropriate correction to use. 
iBhattacharva et all ( 1201 ih find that a correction of Ncor ~ 
N{1 - Ar-°'^5) is more suitable and lLukic etal] ( l2009b note that 
the halo concentration parameter for haloes with an NFW pro- 
file a lso affects the correction, a result that is corroborated by 
iMore et al. (2011)._The more sophis t icated corrections proposed 
bv lEukic et alj ( l2009t) and lMore etal] ( l201lh (the former provide a 
correction for FOF haloes with 6 = 0.2 whereas the latter provide 
a correction valid for different values of b) are based on individ- 
ual halo concentrations. We have not adopted this approach as we 
do not possess profile information for our individual FOF haloes. 
The d ifference between the Warren correction and the lMore et al.l 
ilOl ll) correction (which can be as much as 15% for haloes with 
lower particle counts) becomes slight for haloes with many par- 
ticles. For our smallest, 1000 particle haloes, the correction to the 
Warren formula is of the order of < 5% depending on halo concen- 
trations. Not accounting for this therefore introduces a small sys- 
tematic error which in the very worst case of 1000 par ticle haloes 
would alter the FOF masses by ~ 0.15%. In addition [More et al] 
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Figure 3. The FOF mass function across all simulations and redshifts along 
with our fit (equation [12). The top panel shows the ratio of the data to the 
fitting formula. 



( I2OIII) note that the Warren correction is specifically calculated to 
correct for the bias in the mass function itself whereas their study 
is one that corrects for the masses of individual FOF haloes when 
considered in isolation. 



4.3.3 Sub-Haloes 

As noted in §[3] the AHF halofinder produces halo catalogues that 
contain sub-haloes - bound structures inside host haloes. For the 
purposes of constructing a mass function we have to make a choice 
as to how to deal with these haloes. The CPMSO and FOF halofind- 
ers both identify host haloes only, therefore for compatibility we 
make the natural choice of excluding the sub-haloes from the mass 
functions based on the AHF catalogue. The host haloes detected 
by AHF have properties, including mass, that are calculated based 
on all the bound matter contained within the halo, including the 
sub-haloes themselves. This results in the AHF host haloes being 
directly comparable to the CPMSO ones. 



4.3.4 Low-End Particle Cutojf 

We take a minimuin particle cut-off of 1000 particles per halo in 
this study. This is conservative as often halo-finders report haloes 
with as little as 20 particles, or even fewer. JLiikic et al. (2007) rec- 
ommend a 300 partic le minimum for con structing a mass function 
with a pure PM code. lWarren et al] ( l2006h take an 'aggr essive' 400 
particle minimum in their study and lTinker et alj (|2008|) a 'co nser- 
vative' 400. It has been observed recently bv Reed et al .1 ( 1201^ that 
a minimum particle cut-off of at least 1000 is appropriate. Although 
this latter figure was quoted with tree-based N-body codes in mind 
- rather than a P^M code - we use it to motivate our cut-off of 
1000 particles. We note that for studies in other areas it might be 
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Figure 4. The AHF mass function across all simulations and redshifts. 



acceptable to use haloes with lower particle counts, e.g. in cluster- 
ing studies we need to reliably identify and locate a halo, but not 
necessarily know its mass or other properties precisely. For mass 
function calculations it is important that haloes are assigned masses 
that are as close to their correct values as possible, and more parti- 
cles are required to ensure this is the case. 



4.4 The Universal FOF Function 

We find that the halo mass function based on our FOF halo cata- 
logues follows an universal fit applicable for data from all simula- 
tions across all redshifts. The fit is good to within ~ 10% for most 
data points (Figure [3} and takes following the form: 



^1 +1 



-7/0 



(12) 



where A = 0.282, a = 1.406, /3 = 2.163, 7 = 1.210. This 
fit is valid in the range: —0.55 ^ Incr^^ < 1.31, which at 
2 = corresponds to haloes with masses between 1.8 x 10^^ and 
7.0 X IO^^/i^^Mq. Our largest halo, found in the 3.2 /i^^Gpc box 
at 2 = has a mass of 6.4 x 10^^/i"^Mo and Incr"^ = 1.17, 
whereas our highest Ina"^ value comes from a z = 3 halo in the 
same simulation, with a mass of 1.4 X lO^^/i-^Mg. It should be 
noted that this extreme halo is most likely the result of the FOF al- 
gorithm linking together two large haloes via a bridge of particles - 
a systematic effect known as overlinking, discussed in §[5]below. At 
higher Incr"^ values, i.e. for high mass/high redshift, rare haloes, 
the scatter about the fit increases dramatically due to shot noise. 

Whilst it is a remarkable result that over such a large range of 
z and cr we observe a halo mass function that approximately con- 
forms to a universal shape we urge caution when using this fit for 
two reasons. Firstly the mass function is not completely universal 
and exhibits a modest redshift evolution (see § |5](. Secondly it is 
difficult for a mass function derived from FO F haloes, that are ar- 
bitrar ily shaped, to be related to observables dKravtsov & BorgMiil 
|20l3). 
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Fig ure 5. Evolution of the parameters in equation [12] for CPMSO haloes. 
The lXinker eTal] fcOOSh values ai'e shown in blue for comparison, the model 
used in this work is shown in red. The parameters have been fitted in the 
order: 7 — > A — > o — >■ /9. 



4.5 Mass Function for Spherical Overdensity Haloes 

4. 5. 1 Redshift Evolution 

We show in Figure |4] the mass function from our AHF haloes. In 
contrast to the FOF mass function it is clear that a universal fit is not 
appropriate for spherical overdensity haloes (a fact also born out in 
the results from the CPMSO halofinder. For a comparison between 
the CPMSO and AHF halo mass functions see Figure [7] and dis- 
cussion below). We again adopt the fitting function in equation [T2I 
but with a parameterisation that includes a redshift-dependence, as 
discussed below. Since running AHF on multiple checkpoints from 
our largest simulations is computationally expensive we utilise the 
data from the CPMSO on-the-fly halofinder in order to investigate 
the redshift-dependence of the parameters for the spherical over- 
density case. In ^ l4.5.2l below we discuss the differences between 
the three halo finding methods and re-scale the fit to match the 
AHF results. We perform a least-squares fit on the CPMSO data 
for 36 output reds hifts, initially allowing all of the parameters in 
equation[T^to vary. [Tinker et alj ( l2008l) previously found, based on 
lower-redshift data (0 < z < 2), that the parameter controlling the 
exponential cut-off scale, 7, is approximately constant across their 
range of redshifts. We similarly find that 7 is approximate l y con- 
stant, albeit at a slightly different value from [Tinker et al.l ( [200 8l) 
(see Figure[5] bottom-right panel). Given this result we proceed by 
fixing this parameter at its value in the approximation of universal- 
ity: i.e. 7 = 1.318. We then fit the CPMSO data across all redshifts 
again. We find that using a parameterisation that includes flm.{z), 
the matter content of the universe at a given redshift, enables us to 
capture late time behaviour that differs from that at high-redshift. 
This gives us a fit for A{z) in the following form: 

A{z) = n,^{z) {0.990 X (1 z)-^-'^^^ + 0.074} (13) 

With 7 and A modelled we then repeat the procedure for a and 
then /3. We find for a{z): 
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Figure 6. Ratios between the CPMSO redshift-dependent mass function and the data from the CPMSO halofinder across 36 redshifts from z = — 30. The 
colour scheme is as per Figure[3] 



a{z) = n^{z) {3.136 X (1 + z)-3 058 + 2.349} (14) 
and for P{z): 

I3{z) = n^{z) {5.907 X (1 + z)'^'^^^ + 2.344} (15) 

We show how our model compares to the parameter fitting data in 
Figure [5] 

In Figure |6] we compare our redshift-dependent fit to our 
CPMSO data across all redshifts from z = and z — 30. We 
find that the fit is excellent at low redshifts and remains quite good, 
within 20%, all the way to z ~ 20. For z > 20 our fit still gives a 
reasonable match, but the data has large error bars due to the scar- 
sity of haloes then. Nonetheless, we note that around z — 3 — 4 our 



fit slightly under-predicts the abundances of SO haloes across the 
lower-mass Ina^^ range we cover. There is also an apparent under- 
prediction for lower Ina"^ around z — 8 and there is perhaps an 
over-prediction at very high redshifts (z > 15). 

4.5.2 AHF-based Fits 

Our CPMSO on-the-fly halo finder is by its nature simplified and 
potentially more approximate than AHF. In Figure [7] we show the 
ratio between the mass functions derived using these two SO-based 
halofinders, AHF and CPMSO. In Figure [8] we show the ratio be- 
tween the FOF mass function and the AHF mass function. Both 
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Figure 8. Comparison of the mass functions from the FOF and AHF 
halofinders across redshifts from z = — 26. The colour scheme matches 
the simulations as per Figure[3] 



AHF and FOF exhibit consistent behaviour across all simulation 
volumes and resolutions and for all redshifts. The CPMSO halo 
finder largely agrees fairly well with AHF, typically within 10-20%, 
and much less at low redshifts. There are some systematic differ- 
ences around z = 8 in the transition between the 20 /i~^Mpc and 
114 /i~^Mpc boxes. On the other hand, while both FOF and AHF 
are consistent across all box sizes, there is a systematic trend for the 
FOF to yield more rare, massive haloes than AHF for all volumes 
and re dshif ts. A similar t rend w as noted previously bv lReed et al.l 
l l2007l) and [Tinker et"aL I l l2008l) . Based on these results, we con- 
clude that the AHF (and FOF) data is more consistent across dif- 
ferent box sizes than the CPMSO halofinder. 

Given these differences between the two SO halo finders, in 
addition to the CPMSO fit above, we also provide AHF fits to the 
data, with three different types: 1) a redshift-dependent one, based 
around the redshift evolution we have observed in the CPMSO 
mass function, but re-normalized to the AHF data; 2) a present-day 
one, based on our results at z = 0, for direct comparison with pre- 
vious works; and 3) a more precise high-redshift, "Epoch of Reion- 
isation (EoR)" halo mass function, based on all the AHF data past 

2 = 6. 

Given the relationship between the two SO halofinders in Fig- 
ure [7] we expect that the redshift parameterisation presented in 
equations 11 3ti 1 5 1 applied to AHF data leads to a mass function that 
is roughly correct. The fit to AHF at z = can be improved by 
introducing a slightly different parameterisation for A, as follows: 

A{z) = Q.m{z) {1.097 X (1 + + 0.074} (16) 

The combination of equations 1 141 1 151 [T6l and 7 = 1.318 re- 
sults in the ratios shown in Figure |9] We see that this fit is accurate 
to ~ 10% for redshifts less than z — 15. The shape of the function 
in the low Incr"^ range is slightly different for z > 8, although 



the amplitude and the shape are both correct for Incr ^ > 0.3 at 
2 = 8 - 10. 

We now provide more accurate parameterisations for our AHF 
data for several different redshift ranges. At z — we have: 
A = 0.194, Q = 2.267, P = 1.805 and 7 = 1.287 valid in 
the range —0.55 ^ Incr^^ < 1.05. As the redshift evolution 
of the spherical overdensity mass function is mainly apparent at 
later times (2 < 3) we can provide an universal fit that is appro- 
priate for high redshift studies, for example for probing the EoR 
and the Cosmic Dark Ages. An AHF mass function fitted to all 
our data from z = 6 upwards has the following parameterisation: 
A = 0.563, Q = 0.874, /3 = 3.810 and 7 = 1.453 valid in the 
range —0.06 ^ Inir"^ < 1.24 (corresponding at redshift z = 8 to 
a mass range of 3.5 x 10® - 6.3 x lO^^/i^^M©). These values are 
notably different to those given by our z = redshift parameterisa- 
tion, due to the lack of constraining data for low values of Incr"^. 
The ratio of the EoR fit versus our data is shown in Figure [TO] 
There is considerable scatter in the high Intj"^ end due to data be- 
ing incorporated from high redshifts when haloes are scarce. For 
the range 0.05 ^ lncr~^ < 0.5 the fit is accurate to within 10%. 

4.6 Varying A 

For practical purposes it is useful to have a simple method for ad- 
justing the spherical overdensity mass function fits derived above to 
account for haloes defined with different overdensities. To this end 
we have run AHF on the l/i~^Gpc box using values for A ranging 
between 100 to 1600, at z = 0, 1 and 3. The data has then been fit- 
ted with A incorporated into the parameterisation. This allows, to 
good precision at lower redshifts, a mass function based on a given 
A to be inferred from our Aits fits. The results from this procedure 
are shown in Figure fTTI 
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The parameterisadon is as follows. We assume that the 
/a=i78 mass function is suitable for describing mass functions 
with different choices for A when it is suitably scaled by a function, 

T{A,a,z): 



(17) 



We find that the following form for F is suitable for describing our 
data: 

r(A, a, z) = C{A) {^y^^^ exp [p (1 - ^) /a^] (18) 
Where: 

C(A) = 0.947exp [0.023 ( 1^ - l)] 
d{z) = -0.456f^m(^) - 0.139 

p = 0.072 ^''^ 
q = 2.130 
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Figure 11. The effect of the choice of the overdensity criterion. A, on the 
mass function. For z = 0.1,3 (top panel to bottom respectively) the ratio of 
mass functions with a variety of choices for A are shown relative to the 
Ai78 function. All data has been calculated from the l/i~^Gpc simulation. 
Fitted curves, described in the text, are shown as dashed lines. 



The redshift-dependence of the fit is based solely on Qrn{z), via 
the d{z) parameter. 



5 SUMMARY AND DISCUSSION 
5.1 Fitted Functions 

We have provided a number of fitted mass functions based on our 
data. These are summarised in Table|2] 



5.2 Comparison to Existing Fits 

In Figure [12] we show our universal FOF mass function at z = 
relative to a number of fits available i n the literature, includin g 
the widely-used Press-Sc hechter (PS ) ('Press & Schechtei^ IT974h . 
Sheth-Tormen (ST) ( Sheth & Tormen 2002) and Jenkins et al. 
(Jenkins et al. 2001) fits, as well as a numbe r of more recent fits 
dWarren et aL ,2006> . iReed et all I2OOI I2OOI ICrocce et al.l I2OIOI : 
ICourtin et alj|201 ll ; lAngulo et alTT2012l) . We note that the fits by 
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Table 2. Mass function fitting parameters. For details on tlie redshift evolu- 
tion of these parameters see sections r4.5.1l and l4.5.2l 
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0.282 


0.287 
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0.194 


0.563 
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1.406 


1.478 


1.478 


2.267 


0.874 
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3.810 
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Figure 12. Ratios between our universal fit for FOF haloes and a number of 
FOF fits from the literature for 2 = 0. 
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Figure 1 3. Ratios, at 2 = 0, between the AHF 2 = fit (solid black line) 
versus the lTinker et alj )2008l) fit (solid blue line), the redshift parameterised 
CPMSO fit (dotted red line) and the AHF precise fit adjusted to A = 200 
using eguation flTl edotted black line). 



redshift-dependent mass function compares to the AHF 2: = fit 
and the ratio between the AHF 2; = fit and the A = 200 adjusted 
AHF z = fit (calculated using equation 117b. The iTinker et al.l 
( I2OO8I) fit is based on data in the range -0.6 < Ino-"^ < 0.4, while 
our fit includes data over the range —0.55 < Ino""^ ^ 1.35 across 
all redshifts, or specifically —0.5 5 < Ing"^ ^ 1.05 at z = Q. 

We show an agreement with [Tinker et alj ( I2OO8I) to ~ 5% for 
haloes in the mass range 10^^ - lO^'^/i'^M©. Past ~ IO^'^/i'^Mq 
we predict a low er collapsed mass fr action, although t his is outside 
the ran ge of the [Tinker et al. I j2008l) fitting data. The [Tinker et al.l 
( l2008l) fit is based on A200 haloes and we show how our A200 
prediction via equation [T7] compares to the underlying precision 
AHF plot. We not e that this is slightly less congruent with the 
iTinkeLet a lj (|2008|) fit than the A 178 prediction, although it is still 
in agreement to within 10% for most of the mass range in question. 



iReed et alj (l2007h andlCrocce et alj feoioh are redshift-dependent 



mtt -dep e 

while iPress & Schechteil ( 11974) : ISheth & Tormenl JioOZ) and 



ICourtinetal.U201lh are parameterised using 5c(z) . Aside from the 
classic, and less precise iPress & Schechteil i l 19741) expression, our 
fit agrees to ~ 10% with all the others for most of the mass range, 
and, for very ma ssive clusters (> 1 0^^fe~^M(;)) it matches very 
closely the recent lAngulo et all ( 1201 2i ) fit based on the Millennium 
series of simulations. 

With the exception of the PS, ST and lReedet ( |2003|) re- 
sults, the largest discrepancy both with respect to our fit and overall 
scatter is for very large, rare haloes, whose statistics in most sim- 
ulations is poor. A reliable result in this range requires very large 
simulation volumes of tens to hundreds of Gpc^, comparable to the 
volume of the observable universe. 

As we discussed above, the only oth er robust SO-based halo 
mass function fit is the one provided by [Tinker et all ( |2008|) . In 
Figure 1 131 we show the ratio between our AHF z = fit and the 
[Tinker et alj j2008l) fit at 2 = 0. We also show how the CPMSO, 



5.3 Universality 

The results presented here are in line with other recent studies of 
the halo mass function that have addressed the question of univer- 
sality, i.e. whether the mass function can be considered invariant 
over all redshifts or cosmological models when expressed in a suit- 
able form. It is beyond the scope of this paper to d iscuss invari- 
ance u n der different c o smolo g ies, for work on this see ljenkins et al.l 
( I2OOII): IWarren et all j2006t). [Tinker et all j2008|) and in particu- 



Courtin et al. l 2011 ) and Ishattacharva et al. 201 ih . The mass 



lar 

function has been found to depend only weakly on cosmology 
when couched in terms of FOF haloes with a fixed linking length. 
More study is required to address the question of how SO halo mass 
functions differ in varying cosmologies. 

We find that the assumption of halo mass function univer- 
sality across redshift is approximately valid under certain condi- 
tions, but is violated in general, in qualitative agreem ent with the 
conclusions of other recent studies. [Tinker et al] ( 12008 ) present an 
SO mass function that includes redshift-dependent parameters for 
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z ^ 2.5. Similarly ICrocce et all ( 1201 Ol) and Ishattacharva et al.1 
1 I20T l|) present F OF mass fun c tions with redshift-dependent param- 
eters for 2^2. iReed et aP (l2007h produce a high-redshift mass 
function that contains a redshift-dependence via the effective slope 
of the power spectrum at the scale of the halo radius. We observe 
here that the FOF halo mass function for fixed linking length is 
close to uni versal for a very larg e redshift range (from 2 = 26 to 
the present). I Jenkins et alJfeOOlf) noted that taking a constant link- 
ing length keeps the mass functio n closer to a universal form, a 
result corroborated by iLukic et al] ( 12007.) who found th at even for 
high redshifts the mass function of I Warren et alj ( l2006h was suit- 
able for FOF haloes, despite it being calibrated on 2 = data. 
However, as we discuss in more detail in § 15.4.11 below, using a 
fixed linking length for both low and high redshifts is inappropriate 
and massive FOF haloes often suffer from significant over-linking 
problems, which results in an over-estimate to the halo mass func- 
tion in that range. While we derived an universal fit for the FOF 
halo mass function for fixed linking length of 0.2, there clearly is 
some, albeit modest, redshift evolution about this fit. We show in 
Figure[T4lthe residuals between our FOF universal fit and our data 
for nine different redshifts. The amplitude of the data relative to 
the fit drops away markedly for lower Incr"^ values, to the extent 
that by 2 ~ 8 and above the data is around 20% lower than the 
prediction around Incr^^ ~ 0.5. 

In contrast, the SO-based halo mass functions are clearly not 
universal over the entire redshift range we consider, which was 
also previously noted bv IXinker et alj ( l2008h based on data over 
a smaller range in redshift and a than our data. Here we presented 
SO halo mass function fits that are applicable across a large range of 
redshifts, by explicitly introducing an Q.m- and redshift-dependent 
parameterisation of the mass function. 



5.4 Haloiinder Comparison 

5.4. 1 Relationship Between SO and FOF Haloes 

Care must be taken when using an FOF mass function, especially 
at high redshift, as overlinking of the larger haloes becomes a fac- 
tor. This occurs when one halo is joined to another by a bridge 
of particles that all lie within the linking length, thus creating ab- 
normally massive haloes in the simulation catalogue that are not 
physically meaningful as compact, virialized structures. The dif- 
ferences between an FOF halo and an SO halo, and the respec- 
tive mass functions derived from haloes found using the two algo- 
rithms, are significant. A typical FOF halo is arbitrarily shaped and 
has been described in the literature as a demarkation of an isoden- 
sity contour in real space, with a value of Afof - the overdensity 
of matter at the halo boundary - that is under some debate. Esti- 
mates for Afof include Afof ~ 2b~^, i.e Afof ~ 250 for 

f II ' 

the standard choice of 6 = 0.2, (Navarro et al. 1997); A fof ~ 

0.48by^, i.e. Afof ~ 60 JEacev & Cole 1994; Surnm ers et al.l 
1 19951 ; lAudit et alJI 19981); Afof ~ 74 JWarren et al.ll2006h ; and 
Apnp ~ 81.62 lMoreetai] ( l201lh . all for b = 0.2. In contrast the 



SO algorithm produces haloes that are by definition spherical and 
that adhere to a strict overdensity criterion, e.g. A — 178 times the 
background matter density. This overdensity refers to the average 
density of matter contained within the halo. Haloes are generally 
very centrally-concentrated, thus the overdensity of matter at the 
boundary of an SO halo is much lower than the mean overden- 
sity and the relation be tween the two is d ependent on the profile 
of the halo in question. iMore et al] ( 1201 ih note that for a singular 
isothermal sphere (SIS) density profile, p(r) oc r~^, an overden- 
sity at a halo's boundary of A^dgc ~ 60 corresponds to an enclosed 
overdensity of A^aio ~ 180. This loosely links the value of b in 
the FOF algorithm to A in the SO algorithm -by assuming that the 
FOF halo is spherical- but the simplified SIS profile does not repro- 
duce simulated dark m atter haloes well. More appropriate profiles 
include the TIS profile ('Shapiro et al."l999'; llliev & ShapirolbOOlh 
and the NFW profile. More et al. (2011) used the latter, sampling 
NFW haloes with random particle realisations and show that their 
value of Afof ~ 81.62 corresponds to a range of Ahaio between 
~ 200 —600 depending on the concentration parameter of the NFW 
haloes. iLukic et aT use a similar approach and were able to 

recover, to 5% accuracy, an SO mass function from FOF haloes by 
individually relating the FOF halo masses to their SO counterparts. 
The translation between the two was based on the concentration 
parameter of the haloes' NFW profiles and three empirically fitted 
parameters. 

5.4.2 Choice of Linking Length and Overdensity Criterion 

Given these analyses it is not surprising to observe a marked dif- 
ference in the physical structure of FOF and SO haloes. Figure [Tsl 
shows a large halo from our 20 /i~^Mpc simulation at 2 = 8 cap- 
tured by the CPMSO and the FOF halo finders, respectively. We 
show the FOF halo found with linking lengths of 0.2, 0.15, and 0.1 
in grey, blue and red respectively. This is the second largest halo 
in our volume and it illustrates particularly well the differences be- 
tween the two algorithms. For this halo at this redshift a linking 
length of 0.2 is far too aggressive and significant overlinking has 
occurred. 

In Figure [T6] we compare the halo mass functions based on 
FOF haloes with linking lengths of 0.2, 0.15, 0.1 and SO haloes 
with overdensity choices of 100, 178, 200 and 800, at 2 = 
from the l/i~^Gpc box (left-hand panel) and at 2 = 8 from the 
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overdensity and FOF with a variable linking length. All ratios are plotted on a base of the A = 178 results. 



114/i~^Mpc box (right-hand panel). The z — Q result shows that, 
for a suitable range of o, an overdensity of 178 is comparable to 
a linking length of 0.2 to within 10%. However, for a 0.2 linking 
lengh there appears to be a trend towards lower overdensities for 
the higher mass haloes. This is likely due to the increasing influ- 
ence of overlinking on the masses of the larger FOF haloes. The 
z = 8 results show that a linking length of 0.2 is not at all com- 
plimentary to an overdensity choice of 178 in the EdS regime of 
structure growth. In fact, we see that an overdensity of 100 is more 
consistent with a linking length of 0.2 at z = 8. The equivalent 
linking length for an overdensity of 178 lies between 0.2 and 0.15 
and much closer to 0.15 than to 0.2, and a linking length of 0.1 
roughly corresponds to overdensity of 800. Given these results we 
expect there to be an evolution in the relationship between over- 
densit y an d linking length. T he empirical results of ICourtin et al.l 
( l201lh and lMoreet.Zl l l201lh both contain a redshift dependence 
(the former via A(z ), and the latter vi a the concentration parame- 
ter, c{z)). Whilst the lMore et al.l ( l201~i1) study looked at haloes up to 
2 = 2.5 further study needs to be undertaken to investigate the suit- 
ability of existing relations between linking length and overdensity 
at higher redshifts. 



5.4.3 SO vs. FOF Mass Functions 

The differences we observed above between haloes found using 
the two approaches filter through to secondary results derived from 
halo catalogues, including the mass function. Figure [8] shows the 
systematic difference between the AHF (A = 178) and the FOF 
(b = 0.2) mass functions over a wide range of redshifts. There is 
a close similarity between the two mass functions at z = 0, espe- 
cially for lower mass haloes, with the collapsed fraction of mass 
in higher-mass haloes greater in the FOF case. Gradually, as we 
move to higher redshifts, we see the amplitude of the SO mass 
function fall to around 80% of the FOF mass function at lower 



masses, with a much more marked decrease at higher masses. Past 
redshifts of around z ~ 6 the high mass tail of the SO mass func- 
tion is around 50% lower than that of the FOF mass function. The 
causes of the difference between the two mass functions can be 
summarised as contributions from (1) the relationship between the 
masses of a given SO halo and its FOF counterpart, which depend 
on the choices of b and A and also on the concentration param- 
eter of the haloes; (2) the amount of overlinking of FOF haloes; 
(3) the relative mass difference in the two halo types that arises 
from pseudo mass evolution, as discussed in § 13.41 and (4) other 
systematic effects including the SO algorithm not correctly inter- 
preting the properties of non-spherical haloes and both algo rithms 
failin g to reliab l y desc ribe merging systems. iLukic et a y ( l2009h 
and iMore et al.l (1201 ll) have addressed the first f actor in detail, 
andlDavis et al.l (ll985h;lMa & Bertschinged ( ll995l) ; ICole & LacevI 
l ll996l) ; iLukic et all ( |2009| )~ have investigated the second and, to 
some extent, the fourth factors. The third factor has not been stud- 
ied in detail and remains a topic for future study. 
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Figure 15. Image of a large halo in the 20 Mpc box at z = 8. The 
circle represents the extent of the Aiyg cutoff used in CPMSO. The z di- 
rection has been projected onto the x-y plane. The CPMSO halo mass is 
3.1 X 10^''h~^MQ and it contains 9.3 million particles. The dots represent 
aggregations of at least 20 particles found in the FOF version of the same 
halo. Grey shows the halo captured with a linking length of 0.2, blue 0.15 
and red 0.1. The masses (particle counts) are 4.8x IO^^/i^^Mq (13.1 mil- 
lion), S.TxIO^O/i-IMq (8.9 million) and 2.1 xIO^'^/i-IMq (5.8 million) 
for b = 0.2, 0.15, 0.1 respectively. 
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